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Abstract 

We introduce an event-driven simulation scheme for a simplified model of dense suspensions under 
flow. In this model, hydrodynamic interactions are neglected. The viscous drag force exerted on 
a particle is proportional to the difference between the fixed velocity of the underlying fluid, and 
the velocity of the particle. Our event-driven approach is based on an exactly-derivable equation of 
motion which depends solely on geometric information. Our method allows for a robust extraction 
of the instantaneous coordination of the particles as well as contact-forces statistics and dynamics, 
under any chosen mode of deformation. It can also be used for generating jammed packings under 
shear or compression. 

1. Introduction 

The viscosity of a suspension of immersed hard particles was computed early on by Einstein at 
low volume fraction ^ [H . At larger volume fractions 0, approaching random close packing 0c where 
jamming occurs ji], yj, the viscosity diverges and the motion of the particles becomes increasingly 
correlated 0, [Ej. These observations suggest that dense flows may be governed by a dynamical 
critical point. Understanding this critical behavior would resolve a geometrical question fundamental 
to various problems of condensed matter: how can particles move while still avoiding each other at 
high density? 

Computer simulations are playing an increasingly important role in the investigation of the re- 
lationship between microscopic structure and rheology of amorphous systems - concentrated sus- 
pensions in particular - and are well suited to investigate this question. There exist a number 
of simulational methods for suspensions which take into account hydrodynamic, interparticle and 
Brownian interactions. These methods are reviewed in jl, 0]. In this work we study a model of sus- 
pensions where hydrodynamic interactions are neglected, which simplifies the dynamics considerably. 
This class of models has received a lot of attention recently 0, Isl, 0, 10, 11, 12, 13, 14, 15 1 as they 
appear to capture, qualitatively at least, the observed presence of a jamming transition where the 
viscosity diverges. 
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The model under focus consists of rigid frictionless particles immersed in a flowing Newtonian 
fluid. The flow of the fluid is characterized by an afline velocity field, and the dynamics is assumed to 
be overdamped. Hydrodynamic interactions are neglected, such that the viscous drag forces acting 
on the particles are proportional to the difference between the particle velocities and the velocity of 
the underlying flowing fluid. 

A similar model was initially proposed by Durian jo], [lo[ for the description of sheared foams, 
where soft repulsive interactions are introduced between the constituent bubbles. At packing fractions 
below jamming, the quasi-static limit (i.e. vanishing strain rate) of such a model corresponds to rigid- 
particle dynamics, as the soft particles always have enough time to rearrange to a state of vanishing 
energy, so that particles do not overlap. The quasi-static limit of this model was recently studied 
numerically by imposing vanishingly small stresses 0] or strain rates js], [ll|, 12, 13|. Two methods 

have been employed to simulate this system: (i) small imposed strain increments are followed by 

the minimization of the potential energy 12|, |l3|, or (ii) integration of the overdamped equations of 



motion, where the damping is proportional to the difference between the particles' velocities and the 
affine flow 0,S[il|. 

One difficulty of method {i) is that the energy landscape of soft particles below jamming displays 
many flat directions - the so-called floppy modes. Various minimization methods in quasi-static 
simulations may result in uncontrolled evolution along these floppy modes. The steepest descent 
minimization method does not face this problem, as in principle no undesired motion occurs along 
floppy modes during minimization. However, the usefuUness of the steepest descent method is ham- 
pered by its extremely long running time in the quasi static limit. The integration of the overdamped 
equations of motion |8|, of method (ii) above requires ever-smaller strain rates to allow for enough 
time for the system to reorganize and relax towards the zero-energy states, in order to probe correctly 
the rigid particle limit. In practice the integration step is chosen to be rather large, allowing to probe 
large systems at the expense of loosing accuracy in the microscopic dynamics. Here we introduce 
a new method that allows us to extract precisely which contacts are present, how contacts appear 
and disappear, and how these events affect the contact forces throughout the sample. Our method 
also enforces the correct motion in the space of floppy modes, as steepest descent does, but it is less 
expensive computationally. 

Our approach is based on the exactly derivable equations of motion in the hard sphere limit. These 
equations are used to build an event-driven simulation. The equations of motion are entirely based 
on geometric information, which allows for the calculation of contact forces between the constituent 
hard particles, and hence their distribution and evolution. The main idea is to evolve the system 
according to the equation of motion, while carefully handling the formation of new contacts and the 
opening of existing contacts between particles. Our method shares some similarities with the method 
of contact dynamics 16|] used in dry granular materials, in which contact forces and velocities are 
resolved iteratively under a set of complementarity relations. 

This paper is organized as follows: Sect. [2] describes the simplified model which our simulation 
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method follows. Sect. E] derives the equation of motion for the presented model, building on its main 
assumptions. Sect. S] illustrates the concept of the instantaneous contact network, which is geometric 
information on which the equation of motion is based. Sect. |5] contains an elaborate description of 
our simulational scheme and ends with a perscription for generating jammed configurations under 
shear or compression. Sect. [6] presents results from our simulations, illustrating the utility of our 
event driven approach. Sect. [7] presents concluding remarks. 



2. Model of the flow of hard particles immersed in a viscous fluid 

Consider a system of rigid frictionless spherical particles immersed in a viscous fluid, at a volume 
fraction (j). The fluid phase is assumed to flow with some affine velocity field denoted as V^, and 
the particles are assumed not to alter the flow of the fluid phase. The fluid viscosity is r/o? and the 
Reynold number is negligibly small. 

These assumptions lead to two consequences: 

(i) The fluid exerts drag forces on the particles that are proportional to the difference between the 
particles' velocity and the velocity of the fluid at the position of the particles. Denoting the 
/c*"^ particle's velocity by and its position by -R^, the force exerted on this particle by the 
fluid is (see Fig. [T]) 



(ii) The dynamics is overdamped, therefore forces on each particle are always balanced. Thus, the 
drag force exerted on the /c'th particle by the fluid is counter-balanced by the contact forces 
exerted by the particles in contact with that particle. The particles are assumed to be strictly 
rigid and frictionless, so the contact forces between the particles are purely repulsive. We 
denote the contact forces between the /c*^ and i^^ particles by fik] here and in the following we 
will conventionally set repulsive forces to be positive, i.e. fik > 0. The force-balance condition 
then reads 

— Fk = fik^ik , (2) 

% in contact with k 

where n^fc is a unit vector pointing from the center of the i^^ particle to the center of the fc*^ 
particle. 

In Fig. [1] both results described above are illustrated, see figure caption for details. 



3. Equation of motion 

From relations (i) and [ii) in the previous section we now derive the exact equation of motion. 
Consider a system of N hard spheres (referred to in the following as particles) in a volume V in 
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Figure 1: Illustration of the two main relations defining the model. Rigid particles are immersed in a viscous fluid 
described by the velocity field (horizontal arrows on the left). The fc"^ particle in the center is moving with the 
velocity Vk, and the difference between the k*'^ particle's velocity and the velocity of the underlying fluid is Vk — V^ ■ 
The drag force on the k^^ particle, Ffc, is proportional and opposite in sign to Vk — v ■ Finally, the overdamped 
dynamics imply that the drag force Fk is balanced by the vector-sum of the contact forces /i, ji and /a. 

d dimensions, such that there are no overlaping particles, and some of the particles are exactly in 
contact: the distance between their centers is equal to the sum of their radii. When the system 
is unjammed, the number of contacts in the system remains smaller than the number of spatial 
degrees of freedom Nd — d. Note that we subtract d translations but not rotations due to the periodic 
boundary conditions. We denote the d dimensional vector of the z'th particle coordinates as its 
time derivative as V^, and define the directional differences Rij = Rj — Ri, the pairwise distances 

Tij = \J Rij ■ Rij , and the normalized directions = Rij /rij. We will refer to the contact network 
as the set of all pairs of particles that are in contact at some instance in time, and the geometric 
information that accompanies the network, namely the directions fiij, and the pairwise distances r^j. 

The rate of change induced to a pairwise distance rij, given some vector of particles' velocities Vk, 
follows: ^ 

= Yl 7^ ' = Y^^ik - Sik) nij ■ Vk = {Vj - Vi) ■ fiij . (3) 

k ^^k k 

Fig. ([2]) illustrates how the velocities Vk translate to the rate of change of pairwise distances rij. 
The above relation consists of a linear transformation of vectors from the space of the particles (of 
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Figure 2: Given the velocities of the particles 14, the rate of change in the pairwise distance is the projection of 
the velocities on the direction Eq. 0. 



dimension Nd) to vectors in the space of contacts (of dimension Nc). We define 



S = 



(4) 



then Eq. (I3l) can be written in bra-ket notation as 



f) =S\V) . 



(5) 



Here and in the following we denote vectors from the space of the particles with upper-case notation 
(e.g. \V) for particle velocities), and vectors from the space of contacts with lower-case notation 
(e.g. |r) for the vector of pairwise distances). 

As the particles are completely rigid, they cannot penetrate each other. In addition, a given set 
of contacts does not change except at discrete points in time (see discussion in Sect. H]). Except for 
at those discrete time points, the velocities must satisfy 



We define the nonafjine velocities as the difference between the total velocities of the particles and 
the affine velocity of the underlying fluid: 



r) = S\V) = . 



(6) 



1^°^) = \V) - . 



(7) 
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Combining with Eq. the condition of no overlaps becomes 



S\V) = S {\V^) + IV'''')) = , (8) 

or 

^l^'^'^) = -S\V^) . (9) 

Since \V^) is perscribed by the model, the right-hand-side (RHS) of the above equation is known. 
However, it cannot be solved for jV^'^'^), as S is generally not square, and hence cannot be inverted. 
This is because the number of contacts Nc (the number of linear equations in Eq. (Q to be solved), 
is smaller than the number of degrees of freedom Nd — d (the number of variables in Eq. (jO])). 

We now use the relations in Eq. ([T]) and Eq. ([2]). We first recall that the drag forces \F) exerted on 
the particles by the fluid are proportional and opposite in sign to the difference between the particles' 
velocities \V) and the affine velocity of the fluid l^'^). Thus the drag forces are proportional and 
opposite in sign to the nonaffine velocities defined in Eq. ([7]), namely 

-\F)=Vo{\V)~\V^))=Vo\Vn , (10) 

which is the same equation as Eq. ([1]), but in bra-ket notation. Without loss of generality we choose 
?7o as our unit of viscosity; the above relation is thus 

- \F) = . (11) 

Due to the overdamped nature of our model the drag forces exerted on the particles by the fluid are 
exactly balanced by the contact forces between the particles. Starting from Eq. ([2]), we have 

— Fk = fiknik = fij {Sjk — Sik)nij = fij . (12) 

i in contact with h all contacts i^j all contacts i^j ^ 

The above relation consists of a linear transformation from contact-space to particle-space; the 
transformation is identical to the one defined in Eq. (jl]), namely However, it operates on a 

vector in contact space, as opposed to Eg . (|3l) where the the transformation operates on a vector in 



particle-space. Thus, Eq. ( IT^ becomes [l7|, [l8| : 

-\F)=S^\f), (13) 

where |/) is the vector of contact forces. Combining Eq. (11 II) and Eq. (fT^ . we obtain the relation 
between the nonaffine velocities and contact forces: 

\Vn=S^\f). (14) 
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Inserting the above relation in Eq. 

S\V'''') =SS^\f) =Af\f) = -S\V^) , (15) 

where we have introduced the TV-matrix, defined as 

A/" = 55^ . (16) 

By definition the TV-matrix is a symmetric semi-positive definite matrix, hence unless it is singular 
(see Subsect. 15.31 for discussion), Eq. ( 1T5]1 can be inverted in terms of the contact forces |/) as 

I/) = -Af-^S\V^) . (17) 

As mentioned above, the velocity field of the fiuid is perscribed by the model, therefore the above 
relation completely determines the contact forces. Now, using Eq. f lT^ the nonaffine velocities are 
obtained as 

|\/-) = 5^1/) = -s^^^-'s\v'') . (18) 

This completes the derivation of our equations of motion; the velocities consist of affine and nonaffine 
contributions: 

\V) = \V^) + IV"^) , (19) 

where the affine term \ V^) is perscribed by the model, and the nonaffine part IV'^^) is determined by 
Eq. ([HD. 




Figure 3: The instantaneous contact network is defined at any instance in time by collecting the set of all pairs of 
particles in contact, and extracting the geometric information of the corresponding network. 
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4. Instantaneous contact network 



The equations of motion derived above hold for a single instance in time. They are entirely 
based on the geometric information encoded in the operator S, which is defined on the existing set 
of contacts at that time. Given a configuration of the system we consider all the pair of particles 
in contact as our contact network (see Fig. [3]), and based on this set we calculate the operator S. 
However, this set of contacts is ever-changing, as particles move apart or collide. Nevertheless, the 
events of 'breaking' a contact between two particles, or the formation of a new contact between two 
colliding particles, are instantaneous. As a result, one can define time intervals in which a specific 
contact network remains intact. These time intervals are exactly the periods between any such event 
of contact-breaking or contact-formation, as illustrated in Fig. |H 



Figure 4: A typical signal of the instantaneous average coordination evolving in time. A contact network remains 
intact during the time intervals between events of contact formation/breaking. 

The idea of an instantaneous contact network is central to our simulation method. In the subse- 
quent explanations of the simulation method, we will often refer to the contact network, and to the 
consequence of adding or removing pairs of particles from it. 

5. Simulational scheme 

We now describe our simulation scheme. The main idea is to evolve the system by integrating 
Eq. ( !T9|) in time increments St, until an occurrence of an event in which the contact network changes. 
These instances are inferred by considering the following three conditions at all times: 

(a) The distance between the centers of particles that are in contact is (up to integration errors to 
be discussed in the following) the sum of their radii. 
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(6) Particles that are not in contact do not overlap. 



(c) The contact forces |/) (given by Eq. (fT7|) ) between particles in contact are positive. 

The main simulation scheme is as follows: 

(i) Given a contact network at some instance t in which all contact forces are positive, an in- 
tegration step 6t, and the particles' velocities \V), calculate the time tcoi at which the next 
collision will occur. The calculation of the next collision time is identical to the well-known 
calculation of the next collision time in event-driven hard sphere simulations, see for example 



19|. However, the set of pairs of particles considered for a possible future collision are only 



those which are not already in contact. 

(u) Compare tcoi with t + 6t; then 

(ii.a) If tco\ < t + 6t, integrate Eq. (fT9!) up to tcoi, and add the colliding pair of particles to the 
contact network. This addition may render some of the other contact forces of the network 
negative. If indeed the collision results in negative contact forces, they are removed at 
this point by the procedure described in Subsect. 15.11 

{ii.b) If tcoi > t + 6t, integrate Eq. f lT^ up to t + 6t, calculate the contact forces at the new state 
via Eq. f lT7|) . and check whether any of the contact forces has become negative. Remove 
the contacts with negative forces from the contact network. 

(iii) Step (u) produces an updated contact network with strictly positive contact forces and no 
overlapping particles. The new velocities \V) are then calculated via Eqs. ( |T8l) and (|T9i) . Goto 
(.). 

The integration step 6t in our simulations is set by requiring that the number of integration 
steps between events (new contact formations or vanishing contact forces) is on average larger than 
10. Satisfying this criterion requires gradually smaller integration steps for larger or denser systems. 
Moreover, in our simulations we bound the integration step by requiring St < 10~^. 

In the following subsections we will explain how contact forces, which are rendered negative by 
the formation of a new contact, are filtered out of the contact network, and estimate the errors of 
the simulation. 



5.1. Filtering out contacts with negative forces from the contact network 

As mentioned above, the creation of a new contact may lead to breaking of existing contacts. 
The determination of the new contact network after detection of negative contact forces is carried 
out by repeating two procedures: 
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(i) First, we successively remove the contact with the most negative force from the contact network, 
and re-calculate the contact forces via Eq. fll7p after each such removal. The identity of all 
the contacts removed from the contact network during this procedure is stored. This removal 
procedure is carried out iteratively until the contact network is free of negative forces. Then, 
the total velocities are calculated using Eqs. f|T8|) and f|T9|) . 

{a) With the velocities \V) calculated at the end of procedure (i), we now consider the list of 
removed contacts that we have constructed during procedure (i). For each removed contact, 
we calculate the pairwise velocity 



where the pair {i, j) no longer belongs to the contact network. Note that since we are considering 
pairs of particles which are no longer in contact, constraint ([6]) does not apply, and generally 
Vij 7^ 0. For some of the pairs considered one may find that Vij < which means that the 
particles which were a part of the contact network are now approaching each other after 
breaking the contact between them. This scenario violates the no-overlapping condition, since 
this pair is (idealy) just in contact (it was just removed from the contact network), hence if this 
pair has Vij < 0, it will instantly create an overlap. We treat this scenario explicitly - the pair 
with the most negative Vij is then re-connected and added back to the contact network, and the 
contact forces are then re-calculated. If there remain any negative contact forces {fij < 0), or if 
pairs of particles that were removed from the contact network are still approaching each other 
{vij < 0), return to (i). Otherwise, the final contact network has been found. The uniqueness 
of the contact network obtained using the above procedure is demonstrated in the Appendix 



5.2. Maintaining correctness and error estimations 

The integration errors in distances between pairs of particles that belong to the instantanous 
contact network are of order This can be shown by returning to Eq. (jS]), and writing it for a 
pair i,j as 



thus the velocity differences are orthogonal to the unit vectors riij of pairs of particles in contact, see 
Fig. [5l Although the errors in the pairwise distances are very small, they can accumulate over many 
integration steps. To prevent the accumulation of large errors in cases where the lifetime of contacts 
is long, we utilize the following correction procedure periodically: 

(i) For each pair i,j in the contact network, calculate the pairwise errors 5rij = rij — {di + dj), 
where di is the radius of the i^^ particle. 




(20) 



A. 




(21) 
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Figure 5: The continuous (dashed) particles are depicted before (after) a displacement of order St. The change in 
the distance between two particles in contact is of order r^^ — ^ St^ due to the orthogonality of the total velocity 
differences with respect to the unit vector Eq. (EI]). 

[a) Denoting the pairwise errors as \Sr), we calculate the correction displacement field 

\6R) = -S//\6r) . (22) 

(Hi) At this stage the displacement field \6R) is treated as a velocity field; we then evolve the system 
in the direction determined by \SR), either until the entire displacement is achieved, or up to 
a collision induced by the evolution. Given displacements \SR) we calculate the next ballistic 
collision time tcoi of pairs of particles that do not belong to the contact network, and record 
their identities. 

{iii.a) If tcoi < 1, update the coordinates via |-R) ^ \R) + tco\\5R), and add the colhding pair of 
particles to the contact network. 

{iii.h) If tcoi > 1, update the coordinates via \R) \R) + \5R). 

[iv) If a new contact was formed in step {iii), repeat steps (z) — [in) untill no new contacts are 
formed. 

(f ) Calculate the contact forces, and filter out negative contact forces using the procedure described 
in Subsect. 15.11 
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This correction procedure eliminates the accumulated errors in the pairwise distance between particles 
which are part of the contact network; the correction displacements \SR) can be thought of as the 
response of the system to an inhomogeneous expansion of the pairwise distances of the contacts 
within the contact network. The frequency of application of the correction procedure can be tuned 
to achieve the desired accuracy of the pairwise distances within the contact network. 

5.3. Generating jammed configurations 

Our simulation can be used to generate jammed configurations, under shear deformation or under 
compression (or both). Before explaining how such configurations can be generated, we first clarify 
the distinction between jammed and potentially jammed configurations. 

In a potentially jammed configuration the ratio between the number of contacts and the number 
of degrees of freedom is such that the matrix S must have a kernel of dimension one. Thus there 
exists a non-trivial vector of contact forces |/o) 7^ that satisfies force balance, i.e. that solves the 
equation S^\fo) = 0. However, a potentially jammed configuration may not be jammed, since the 
vector 1/0) that solves iS^|/o) = may have negative components, which cannot represent contact 
forces between rigid, frictionless particles. Therefore, in a jammed configuration there must exist a 
vector of positive contact forces |/o) > 0, that solves the equation S^\fo) = 0. 

In order to generate such jammed configurations, one must handle the problematic cases in which 
the system becomes potentially jammed. In these cases, the equation S'^\f) = has a non-trivial 
solution, which means that TV = SS'^ has a zero eigenvalue, which in turn implies that Eq. ( |T7I) 
cannot be solved, thus the contact forces and the velocities cannot be calculated. 

The detection of potentially jammed or jammed configurations requires the tracking of the co- 
ordination number of the constrained subset of the system. The constrained subset is the set of 
particles which are in contact with at least d+1 other particles from the constrained subset (where 
d is the spatial dimensionality of the system considered). The calculation of the constrained subset 
must be performed iteratively, since the removal of some of the particles from the constrained subset 
may exclude other particles connected to them from the constrained subset as well. From a physical 
perspective, the monitoring of the constrained subset elliminates the problems of counting degrees of 
freedom associated with rattlers or particles which are just loosely connected in the contact network. 

With the correct detection of the constrained subset in hand, we next define the constrained 
coordination number z as 

z^i (2(^-1)+ Yl ^0 ' (23) 

\ 2 G constrained subset / 

where Zi is the number of contacts between the i^^ particle and other particles j which have Zj > d+1, 
and N is the number of particles i with Zi > d + 1. The system is potentially jammed if ^ = 2d. 
Eq. ( l23l) is obtained by comparing the number of contacts in the constrained subset | Yli 
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number of degrees of freedom Nd — d+1. The addition of one to the number of degrees of freedom 
assures that S has a kernel of dimension one if the system is potentially jammed. Definition 
should be calculated using integer arithmetics, since it is an exact indicator of isostaticitjfl. 

Once a potentially jammed configuration has been detected, the vector |/o) which solves 5^|/o) 



is calculated [20[. If all its components are positive, the system is jammed. Otherwise, the procedure 
described in Subsect. IS.ll for filtering negative contact forces out of the contact network is employed, 
assigning the solution |/o) as the contact forces. Within the procedure, the constrained coordination 
number must always be tracked; whenever the system is detected to be potentially jammed, the 
solution l/o) should be re-calculated and considered instead of the contact forces. The procedure ends 
with either (z) some contacts being removed such that the system is no longer potentially jammed, 
then the forces and velocities can be calculated as usual, or (u) a contact network being formed with 
a solution |/o) to S'^lfo) = 0, having only positive components, then a jammed configuration has 
been found. 
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Figure 6: Segment of a typical stress vs. strain signal from our event-driven simulation, at volume fraction of = 0.83. 
The stress relaxes during the time intervals in which the contact network remains intact. The sudden upward jumps 
in stress signal the collision of particles. Our simulation method enables an accurate study of the statistics of the 
collisions and their effect on the mechanical response of the system. 



^This indicator is exact if there are no over-constraints regions in the sample. This is the case for polydisperse 
particles where crystalline over-coordinated regions cannot appear. 
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6. Results 



We simulate systems of discs in two dimensions, under simple shear flow. Our system consists of a 
binary mixture, where half of the discs are small and half are large. We add a uniform polydispersity 
of 3% to avoid hexagonal patches of identical particle^. Our unit of length A is the mean diameter 
of the small particles. The large particle have a mean radius of 0.7A and the small particles have a 
mean radius of 0.5A. The velocity field of the underlying fluid is V^{R) = yx. Simple shear flow is 



imposed homogeneously throught the simulation cell using Lees-Edwards boundary conditions [19 
Detection of newly formed contacts is acheived by maintaining data structures of neighbor-lists which 
are updated using cell-subdivisions that allow for the calculation of the neighbor-lists in a time which 
is linear in the system size [19|]. The equation of motion f[T^ is integrated by solving Eq. f lT7|) using 
a standard conjugate gradient algorithm [211]. 

In Fig. we plot a small segment of a typical stress vs. strain signal from our simulations at a 
volume fraction of = 0.83. The stress a is measured using information of the contact forces via 

.^-mn, (24) 

where V is the volume of the system. The signal displays a series of intervals in which the stress 
relaxes, interrupted by abrupt increases in the stress, which are the consequence of collisions (new 
contacts are formed). 

In Fig. ([7]) we display a realization of our sheared systems at = 0.83. The thickness of the 
lines connecting the centers of neighboring particles are proportional to the contact force between 
the particles. 

6.1. Comparison with steepest descent simulations 

An alternative to integrating our equation of motion, Eq. flTIJl) . one can apply small strain-steps 
followed by the minimization of a soft, repulsive interaction potential by means of steepest descent, 
i.e. according to the first order equation of motion 

Rk = -Vo^kU , (25) 
where a possible choice of the potential energy is U = Ylii j>i't'ijij)i and 



[nj) = { dpk-J ' "^iJ < + , (26) 

, ra > di + dj 



^Such patches can be hyperstatic, i.e. over-constrained. Physically these patches move a solid blocks, until they 
break. In our approach these patches are difficult to deal with (although it can be done), as they lead to zero 
modes in the spectrum of A/", which renders Eq. (|17p difficult to solve numerically. We avoid this problem by adding 
polydispersity. 
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Figure 7: A snapshot from our two-dimensional simulations of bi-disperse discs under shear flow. The thickness of the 
lines connecting the centers of the particles represent the magnitude of the contact forces between the particles. The 
colors distinguish between the two sizes of discs. 

where e is an energy scale and di is the radius of the i^^ particle. The equivalence between the 
two approaches stems from (i) the identical dissipation mechanism, which only acts along the the 
nonaffine part of the motion, and {ii) the total evolution that keeps the potential energy at zero since 
the system is studied in the floppy regime, where it is always possible to re-arrange the particles to 
a state of zero potential energy. This is analogous in our simulation to the no-overlap condition 
which is always satisfied. To validate our simulational scheme, we compare the results of our scheme 
with the steepest descent (SD) simulations. In the SD simulation, we construct a tentative contact 
network by considering pairs of particles that satisfy = di + dj + 6r with 6r = 10~^. This network 
is constructed only for the contact force calculation from which the stress is deduced. We choose 
the strain increment for the SD simulations to be ^7 = 10~^ and choose the stopping condition 
for the minimization to be max (Vjf/) < 10~^£:/A. The integration step for the steepest descent 
minimizations was set to be 6t = 0.1. In Fig. [H] we plot the trajectories of the stress vs. strain 
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for both simulation schemes, for systems of = 100 particles at the volume fraction = 0.82. 
The agreement between the two signals vindicates our event-driven scheme. We note that for strain 
increments of ^7 = 10~^, with which the events of contact formations can be clearly singled out 
in the SD simulations, the running time is extremely long, hence the choice of small systems for 
comparison with our method. It took over 30 hours on a conventional workstation to produce the 
data of Fig. [S] using the SD method, whereas it takes just a few minutes with our event driven code 
on the same machine. 
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Figure 8: Comparison between stress vs. strain trajectories from our event-driven simulation and from the steepest 
descent method. Starting from the same initial configurations, the trajectories remain essentially identical for tens of 
percents of strain, vindicating our event-driven scheme. 



7. Discussion 

We have introduced an event-driven simulation scheme for a simple model for dense suspension 
flow. The simulation is based on an equation of motion which is exactly derivable from the model 
assumptions. The simulation maintains the conditions of no overlaps and strictly positive contact 
forces between the rigid particles at all times. We have introduced a procedure that allows for 
the correction of the integration errors in the pairwise distances between particles that are part of 
the instantaneous contact network. We have also provided a perscription for generating jammed 
configurations employing our method, under compression or shear. 

We have demonstrated that our method produces essentially identical results to steepest-descent 
minimization methods, though at greatly reduced computational cost as well as increased robustness 
and accuracy of the contact network information. Indeed, the main virtue of our event-driven method 
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is the control over the identity of the contact network that allows for careful study of the statistics of 
contact formation and breakage. In our approach there are no uncertainties nor arbitrary thresholding 
associated with the definition of which particles should be considered as contacts. This advantage 
becomes increasingly important closer to jamming where contacts rearrange very rapidly, and where 
many particles are very close but actually not in contact. 
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Appendix A. Uniqueness of the contact network upon new contact formation 

The algorithm described in Subsect. 15.11 begins with the full set of contacts, including a newly 
formed contact. We denote this full set of contacts by Q. Upon termination, the algorithm produces 
a subset of contacts, denoted here by A, which satisfies the two following conditions: 

(i) all of the contact forces calculated on the subset A of contacts are positive, and 

(ii) all of the pairs of particles which were just removed from the contact network, are moving away 
from each other. 

In this subsection we prove that the contact network found by our algorithm is unique. We first 
re-iterate that the operators S and (Eqs. (jl]) and ( IT6|) respectively) are generally defined on a 
given set of contacts. When a new contact is formed, the set of pairs on which these operators are 
defined grows by one pair - the new contact. However, as mentioned above, the formation of a new 
contact may render other contact forces negative, in which case those pairs with negative contact 
forces are removed from the contact network. Removing contacts from the contact network re-defines 
the operators S and Af, since the space of contacts has changed. Given the subset of contacts ACQ, 
we denote the corresponding operators defined on this set as Sa and Ma- As opposed to the space 
of contacts, the space of particles never changes; one can, in principle, calculate the full vector of 
velocities given any subset of the full contact space. For instance, the velocities calculated using the 
subset A are 

\Va) = \V^) + \Vr) = 1^")- SIMM'S a\V^) . (A.l) 

We now present the general relation between contact forces and velocities; starting from Eq. f|T^ . 
we consider the full contact network, then 

S^\f) = IV"^) = \V) - \V^) . (A.2) 
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Operating with S on both sides of the above equation, 

SS^\f)=Af\f) = S\V)-S\V'') . (A.3) 

Now, inverting to solve for the contact forces 

I/) = M~^S\V) - Ar-^S\V^) . (A.4) 

Given any velocity vector \V), the resulting contact forces can be obtained via Eq. flA.4p . This 
equation is more general than Eq. f[T7|) . which is obtain from the above relation by demanding that 
S\V) = 0, which is the condition for maintaining an intact contact network. In our discussion we 
relax this condition for a subset of the contact network, since removing contacts will result in non-zero 
pairwise velocities. 

We next note that a contact force which is identically zero has no effect on the dynamics. This 
can be understood by re-examining Eq. fll4p : a zero contact force does not participate in the vector 
sum of the contact forces on each particle, which, in turn, is equal to the nonaffine velocities. This 
allows us to carry out the entire discussion in the full space of contacts Q; the vector of contact forces 
in the full space of contacts, obtained with the velocities calculated on the subset A, is 

\fA)=Ar-'s\VA)-Ar-'s\v'') (a.s) 

Note that the dimension of |/^) is the size of the full set of contacts Labeling contact components 
hj a E Q, and denoting by \6a) a vector which is unity in the a component and zero otherwise, the 
condition of positive contact forces now reads 



{SalfA) > if aeA, 
(<5„|/a)=0 if aen\A. 



(A.6) 



The vector of pairwise velocities induced by the particles' velocities \Va) Eq. (lA.ip is 

\va)=S\Va). (A.7) 

The condition according to which pairs of particles that were just removed from the contact network 
should be moving apart can be written in our framework as 



{Sa\vA) = if aeA , 
{Sa\vA) >0 if a en\A . 



(A.S) 



We now assume that there exists another different subset B (1 B ^ A, and an associated 
velocity vector \Vb) = \V^) — SJ^M^^SbIV^) , the contact forces and the pairwise velocities \vb), 
which also satisfy flA.6p and f lA.Sp for the subset B. We consider the contraction 

{vA-VB\Af-'\vA-VB)>0 , (A.9) 
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which is positive since M is positive definite. Combining \v) = S\V) with Eq. (lA.4p we have 

I/a - fn) = ^f-'S\VA) - ^f-'S\VB) = M-'\va - vb) . (A.IO) 
Inserting this back into Eq. (lA.9p . we find 

{VA - vb\M-'\va - Vb) = {VA - vbUa - Jb) = (vAlfA) + {vbIJb) ' - (vbIJa) ■ (A.ll) 

From the conditions in Eq. (]A.6p and Eq. (lA.Sp . both {vaI/a) = and ("WbI/b) = 0. The inequahty 
(lA.gp now reads 

0<{vA- VB\Af-'\vA - Vb) = -{vAlfs) - (vbIIa) < , (A. 12) 

which is a contradition, since the cross terms {vaI/b) > and (vbI/a) > 0, as all components of 
va,vb, Ja, Jb are non-negative by construction. Hence, there cannot exist another distinct subset 
B A which satisfies Eqs. (1A.6P and (lA.Sp . and therefore the subset A is unique. 
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